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Abstract. We are concerned with central differencing schemes 
for solving scalar hyperbolic conservation laws arising in the 
simulation of multiphase flows in heterogeneous porous me- 
dia. We compare the Kurganov-Tadmor (KT) [3] semi-discrete 
central scheme with the Nessyahu-Tadmor (NT) [57] central 
scheme. The KT scheme uses more precise information about 
the local speeds of propagation together with integration over 
nonuniform control volumes, which contain the Riemann fans. 
These methods can accurately resolve sharp fronts in the fluid 
saturations without introducing spurious oscillations or exces- 
sive numerical diffusion. We first discuss the coupling of these 
methods with velocity fields approximated by mixed finite el- 
ements. Then, numerical simulations are presented for two- 
phase, two-dimensional flow problems in multi-scale heteroge- 
neous petroleum reservoirs. We find the KT scheme to be con- 
siderably less diffusive, particularly in the presence of high per- 
meability flow channels, which lead to strong restrictions on 
the time step selection; however, the KT scheme may produce 
incorrect boundary behavior. 
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1. Introduction 

We are concerned with high resolution central schemes for solving 
scalar hyperbolic conservation laws arising in the simulation of multi- 
phase flows in multidimensional heterogeneous petroleum reservoirs. 

Many of the modern high resolution approximations for nonlinear 
conservation laws employ Godunov's appoach [3S| or REA (recon- 
struct, evolve, average) algorithm, i.e., the approximate solution is 
represented by a piecewise polynomial which is Reconstructed from 
the Evolving cell Averages. The two main classes of Godunov meth- 
ods are upwind and central schemes. 

The Lax-Friedrichs (LxF) scheme [33J is the canonical first order 
central scheme, which is the forerunner of all central differencing 
schemes. It is based on piecewise constant approximate solutions. It 
also enjoys simplicity, i.e., it does not employ Riemann solvers and 
characteristic decomposition. Unfortunately the excessive numerical 
dissipation in the LxF recipe (of order (9((AX) 2 /At)) yields poor 
resolution, which seems to have delayed the development of high 
resolution central schemes when compared with the earlier develop- 
ments of the high resolution upwind methods. Only in 1990 a second 
order generalization to the LxF scheme was introduced by Nessyahu 
and Tadmor (NT) [27]. They used a staggered form of the LxF 
scheme and replaced the first order piecewise constant solution with 
a van Leer's MUSCL-type piecewise linear second order approxima- 
tion The numerical dissipation in this new central scheme has 
an amplitude of order 0((AX) 4 /At) . When applying these meth- 
ods to multiphase flows in highly heterogeneous petroleum reservoirs 
or aquifers we need to use decreasing time steps as the heterogene- 
ity increases, yielding greater numerical diffusion. Kurganov and 
Tadmor (KT) [3] combined ideas from the construction of the NT 
scheme with Rusanov's method [SB] to obtain the first second order 
central scheme that admits a semi-discrete formulation which is then 
solved with an appropriate ODE solver. The resulting scheme has a 
much smaller numerical diffusion than the NT scheme. Due to the 
semi-discrete formulation, this numerical diffusion is independent of 
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the time step used to evolve the ordinary differential equation. This 
property guarantees that no extra numerical diffusion will be added 
if the time step is forced to decrease. For this reason, the application 
of this central scheme results in a new numerical approach to model 
two-phase flows with a much lower numerical diffusion, even in the 
presence of a highly heterogeneous porous media. 

The goals of this paper are (i) to discuss the coupling of NT 
and KT schemes to velocity fields approximated by Raviart-Thomas 
mixed finite element method (See [30J), and (ii) to compare the KT 
semi-discrete central scheme with the NT central scheme for numer- 
ical simulations of two-phase, incompressible, two-dimensional flows 
in heterogeneous formations. Both methods can accurately resolve 
sharp fronts in the fluid saturations without introducing spurious 
oscillations or excessive numerical diffusion. 

Our numerical experiments indicate that the KT scheme is consid- 
erably less diffusive, particularly in the presence of viscous fingers, 
which lead to strong restrictions on the time step selection. On the 
other hand the KT scheme may produce incorrect boundary behav- 
ior in a typical two-dimensional geometry used in the study of porous 
media flows: the quarter of a five spot. 

Numerous methods have been introduced to solve two-phase flow 
problems in porous media. Among eulerian-lagrangian procedures 
we mention the Modified Method of Characteristics [251 [29], the 
Modified Method of Characteristics with Adjusted Advection [22], 
the Locally Conservative Eulerian Lagrangian Method [23] and Euler- 
ian Lagragian Localized Adjoint Methods [26]. Additional tech- 
niques, to name just a few, include higher-order Godunov schemes 
[TO] , the front-tracking method [7], the streamline method [321 G2] 
the streamline upwind Petrov-Galerkin method (SUPG) [H H] and 
a second-order TVD-type finite volume scheme [31] (this procedure 
aims at the modeling of flow through geometrically complex geo- 
logical reservoirs). Each of these procedures has advantages and 
disadvantages. We refer the reader to [281 123] and references cited 
there for a discussion of these methods. 
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We remark that central schemes are particularly interesting for the 
numerical simulation of multiphase flow problems in porous media 
because they have been formulated to solve hyperbolic systems; this 
is not the case for several of the procedures mentioned above, which 
have been developed only for scalar equations. 

Moreover these central schemes were also used to deal with many 
other applied problems: to solve Hamilton- Jacobi Equations (see 
[TT] and [2]), to model the two-dimensional magnetohydrodynam- 
ics (MHD) equations and to study the Orszag-Tang vortex system, 
which describes the transition to supersonic turbulence for the equa- 
tions of MHD in two space dimensions (see [9] and [2D]), to mention 
just a couple of them. 

This paper is organized as follows. In Section [21 we discuss our 
strategy for solving numerically the model for two-phase, immisci- 
ble and incompressible displacement in heterogeneous porous media 
considered here. In Section (3[ we discuss the application of central 
differencing schemes to porous media flows. In Section HI we present 
the computational solutions for the model problem considered here 
and our conclusions. 



2. Numerical approximation of two-phase flows 

We consider a model for two-phase immiscible and incompressible 
displacement in heterogeneous porous media. The governing equa- 
tions are strongly nonlinear and lead to shock formation, and with or 
without diffusive terms they are of practical importance in petroleum 
engineering [T5], [21] • See also [2lj and the references therein for re- 
cent studies for the scale-up problem for such equations. 

The conventional theoretical description of two-phase flow in a 
porous medium, in the limit of vanishing capillary pressure, is via 
Darcy's law coupled to the Buckley-Leverett equation. The two 
phases will be referred to as water and oil, and indicated by the 
subscripts w and o, respectively. Without sources or sinks and ne- 
glecting the effects of capillarity and gravity, these equations read 
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(See pj)] for more details) 



(1) 



V- v = 



v = — 



\(s)K(x)Vp, 



(2) 



ds 
di 



+ V-(/(s)v) = 0, 



Here, v is the total seepage velocity, s is the water saturation, K(x) 
is the absolute permeability, and p is the pressure. The constant 
porosity has been scaled out by a change of the time variable. The 
total mobility, A(s), and water fractional flow function, f(s), are 
defined in terms of the relative permeabilities k ri (s) and phase vis- 
cosities Hi by 



2.1. Operator splitting for two-phase flow. An operator split- 
ting technique is employed for the computational solution of the 
saturation equation (j2J) and the pressure equation ([I]) in which they 
are solved sequentially with possibly distinct time steps. This split- 
ting scheme has proved to be computationally efficient in producing 
accurate numerical solutions for two-phase flows. We refer the reader 
to (22] and references therein for more details on the operator split- 
ting technique; see also [HI [HJ [T7J [18] and [5] for applications of this 
strategy to three phase flows taking into account capillary pressure 
(diffusive effects). 

Typically, for computational efficiency larger time steps are used 
for the pressure- velocity calculation (Equation [1]) than for the con- 
vection calculation (Equation [2]). Thus, we introduce two time steps: 
At c for the solution of the hyperbolic problem for convection, and 
At p for the pressure- velocity calculation so that At p > At c . We 
remark that in practice variable time steps are always useful, espe- 
cially for the convection micro-steps subject dynamically to a CFL 
condition. 

For the pressure solution we use a (locally conservative) hybridized 
mixed finite element discretization equivalent to cell-centered finite 
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differences [ZTJ |22] , which effectively treats the rapidly changing per- 
meabilities that arise from stochastic geology and produces accurate 
velocity fields. The pressure and Darcy velocity are approximated 
at times t m = mAt p , m = 0, 1, 2, ... . The linear system resulting 
from the discretized equations is solved by a preconditioned conju- 
gate gradient procedure (PCG) (See [22] and the references therein). 
The saturation equation is approximated at times t™ = t m + nAt c 
for t m < t™ < t m+1 . We remark that we must specify the water 
saturation at t — 0. 

3. Central differencing schemes for porous media 



In this section, we shall study the family of high resolution, non- 
oscillatory, conservative central differencing schemes introduced by 
Nessyahu and Tadmor (NT) and Kurganov and Tadmor (KT). They 
will be applied to the numerical approximation of the scalar hyper- 
bolic conservation law modeling the convective transport of fluid 
phases in two-phase flows. For the associated elliptic problem (Eq. 
(CE])), we use the lowest order Raviart-Thomas [30J locally conserva- 
tive mixed finite elements. These central schemes enjoy the main 
advantage of Godunov-type central schemes: simplicity, i.e., they 
employ neither characteristic decomposition nor approximate Rie- 
mann solvers. This makes them universal methods that can be ap- 
plied to a wide variety of physical problems, including hyperbolic 
systems. In the following sections we will discuss the main ideas 
of the NT and KT central schemes coupled to the mixed finite el- 
ement discretization mentioned above. We will not repeat here all 
the details involved in the development of the NT and KT schemes; 
instead, we refer the reader to [27] and [3] for this material. 

3.1. The Nessyahu- Tadmor scheme for two-phase flows. Con- 
sider the following scalar hyperbolic conservation law, 
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(3) 



ds 
di 
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subject to prescribed initial data, s(x, y, 0) = So(x,y). Here "v = 
^(x, y,t) and % = y,t) denote the x— and y— components of 
the velocity field v (see Eq. CD). To approximate (j3J) by the NT 
scheme, we begin with a piecewise constant solution of the form 

J2j,k S lkXj,k( x ^y)^ where Sj,k '■= s ( x jiVkit7) is the approximate 
cell average at t — t™ associated with the cell Cj^ = Ij x Ik — 
[xj-i/ 2 ,Xj + i/ 2 ] x [yk-i/2,yk+i/2] and y) is a characteristic func- 

tion of the cell Cj^- 

We first reconstruct a piecewise linear approximation of the form 

s(x,y,t™) = ^2S^ k (x,y)xjA x ^y) 

j,k 



(4) 



j,k 



°j,k 



AX 



(x — Xj) 



Xj,ki x ,v) 



Xj-l/2 <X< X j+1/2 , Vk-1/2 <y< Vk+1/2- 

In Eq. 01]), the discrete slopes along the x and y directions satisfy 
(5a) = jU(:r = Xj , y = y ht t") + O(AX) 

(5b) S^- = ^-s(x = Xj , y = y k , t K ) + O(AY), 

to guarantee second-order accuracy. 

The reconstruction (j3J) retains conservation, i.e.: 

Let {s(x, y,t),t > t™} be the exact solution of the conservation 
law (El), subject to the reconstructed piecewise-linear data (BJ at 
time i = t"\ The evolution step in the NT scheme consists of ap- 
proximating this exact solution at the next time step t — t™ + At c , 
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by its averages over staggered cells, C j+ i/ 2 ,k+i/2 ■— Ij+1/2 x 4+1/2- 
See dashed grid in Figured] (denote n + 1 := t™ + At c ). Let 



S]X\, k+ \ = AXAY / s(x,y,t^ + At c )dxdy 



(7) = -TTFTT7 / s{x,y,t™ + At c )dxdy. 



AXAY 



in, 



Xj < x < x j+ i, yk <y < yk+i 



These new staggered cell averages are obtained by integrating the 
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Figure 1. Evolution step at each time level i™, t m < 
t™ < t m+1 , for the two-dimensional NT central differ- 
encing scheme. 



conservation law (j3J) over the control volumes Cj + i/ 2 ,it+i/2 x KT> + 
At c ] following the same manipulations as described in [19] (denote 
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rv = Ak. anH rv — At °V 
— AX y — Ay ' ' 



At c 



D i+i/2,fc+i/2 - axAY 



s(x, y, t™ + At c ) dxdy 



3+1/2, fc+1/2 



ft, 



AXAF 



AXAF 



t™+At cf y k+1 



Hi, 



x v(x j+1 ,y,r) f(s(x j+1 ,y,r) 



- x v(xj,y,r) f(xj,y,r) 



dy dr 



t™+At c rx j+1 



v v(x,y k+1 ,r) f(s(x,y k +i,T) 



dxdr\ . 



y v(x,y k ,r) f(x,y k ,r) 



The cell average J s(x, y, t™+At c ) dxdy has contributions from 

°j + l/2,fc + l/2 

the four cells C jik , C j+1)k , C j+1)k+1 , and C jjk+1 : 



s(x,y,t™) dxdy 



j + l/2,fc+l/2 



^ fc (x,J/) 



Cj + l/2,fc + l/2 n Ci,fe Cj + l/2,fc+l/2 n Cj ! i ; _). 



(9) 



+ / S« +ljk (x,y) + / S? +1M1 (x,y) 

( <7 + l/2,fc+l/2 n Ci+l,fc Ci+l/2,fc+l/2 n Cj + l,fc+l 



Computing these integrals exactly yields 



S 



( Q K i _i_ C K _i_ Q 1 ^ ^ 



+ 



16 



(10) 



To approximate the four flux integrals on the right hand side of 
(JSj), we use the second-order rectangular quadrature rule for the spa- 
tial integration and the midpoint quadrature rule for second-order 
approximation of the temporal integrals. For instance, letting k+1/2 
be C + At c /2, 



a , 



AXAY 



11a) 



t™+At cf y k+1 

x v{x j+1 , y, t) f(s(x j+1 , y, r))dydr ; 

x l,«+1/2 f f K+l/2x , - K+l/2 f( K+\/2 , 
u j+l,k J ^j+ltk I "i u j'+l,fcH-l ^ 
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AXAY 



y v(x, y k+1 , t) f(s(x, y k+1 , r))dy dr 



a,, 



(Hb) « f f&) + ^Z 2 +1 f(s%& v 

Since these midvalues are computed at the center of the cells, Cj, k , 
where the solution is smooth, provided an appropriate CFL con- 
dition is observed, we can use Taylor expansion together with the 
conservation law ([3]) to get 

(12) ^ = Sj,k ~^ Xvl j,k{fj,k) 7£~ VvK j,k{fj,k) ■ 

Here, (f* >k ) and (fj k ) are one-dimensional discrete slopes in the x 
and y directions, respectively. They satisfy the conditions 

(13a) = ±f( 8 {x = xj, y = y k , t K )) + O(AX) 



(13b) U]&_ = d_ f{s(x = Xjj y = ^ fK)) + G(Ay)) 



in order to produce a second order scheme for the approximation 
of (J3j) . To avoid spurious oscillations, it is essential to reconstruct 
the discrete derivatives given by Equations ([5]) and (|T3|) with built-in 
nonlinear limiters. In this work we use the following MinMod limiter 



( S x)j,k ~ MM ®~K^ { S j-l,k> S j,k> S j+l,k} 

(14a) := MM ( q AS "+W AS "-W ~ AS Uw AS "-i,2,k 
v 1 \ Ax 2Ax Ax 

(14b) := MM (fl^+W 

where A is the centered difference, ASj +1 , 2 k = Sj +1 ^ k — S^ k . We 
refer the reader to [27] and [3] and the references therein for the 
various options for the form of such discrete derivatives. 
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In our sequential scheme, when solving for the saturation in time, 
the total velocity v is given by the solution of the velocity-pressure 
equation. Recall that the solution of Eq. (OQ) is approximated the 
lowest order Raviart-Thomas mixed finite element method. Thus, 
the computed total velocity v is discontinuous at the vertices of the 
original non-staggered grid cells. This constitutes a difficulty for the 
staggered scheme (fSlh which requires the values of the total velocity 
v at these vertices at every other time step. To avoid this difficulty 
we use the non-staggered version of the NT scheme. 

To turn the staggered scheme (IE]) into a non-staggered scheme, 
we re-average the reconstructed values of the underlying staggered 
scheme, thus recovering the cell averages of the central scheme over 
the original non-staggered grid cells. First we reconstruct a piecewise 
bilinear interpolant at the time step n + 1 := t™ + At c 
(15) 



QK+1 | 

°j+l/2,k+l/2 [ 



I qk+1 
\ J j+l/2,k- 



1/2) 



D j+l/2,fc+l/2 -r ^ x 

. (^J + lAfe+l^)/ s 

+ ^ ~ Vk+1 / 2 > 



(x - x j+1/2 ) 



Xj < x < x j+1 , y k <y < y k +u 

as in Equation (jl]), through the staggered cell averages given by (jSJ), 
and re-average it over the original grid cells, giving the following 
non-staggered scheme: 

1 /-p^K,+ l -pjK-f-1 \ 

bj,k = 7 Wj-l/2,fc-l/2 + ^j-l/2,fc+l/2 + ^j+l/2,k-l/2 + ^j+l/2,k+l/2J 



+ 



+ 



16 i y i-v**-y* ' 

L 

/2,/c— 1/2 



I QK+1 



-l/2,k+l/2> 



j+l/2,k+l/2) 



16 L i-iAAs-1/2; 

/ qK+1 

\°j+l/2,k-l/2 



-l/2,fc+l/2/ 

1 qk+1 ~\ 

\ J j+l/2,k+l/2) 



3.2. The Kurganov-Tadmor scheme for two-phase flows. The 

first multidimensional extension of the KT scheme was presented in 
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[3]. This extension used the dimension by dimension approach, that 
is, the numerical fluxes computed along the x and y directions are 
viewed as a generalization of the one-spatial-dimension numerical 
fluxes. This approach consists of the following steps: at each time 
step and at each cell Ij^, 

(i) Compute the difference of the one-dimensional numerical flux 
in one spatial dimension in the x direction keeping y constant 
and equal to y k - Denote this difference by 



(16) 



•^7+1/2, k\y> 



TSfX 

n j+l/2,k 



TJX 

n j-l/2,k 



(t) 



AX 



The numerical flux Hj +1 , 2k (t) is 



Hj+l/2,k(t) 



where 



'v j+ i/2,k{t) f(S+ +1/2 k (t)) 

+W*)/(W(')) 

a j+l/2,fc( t ) 



^•+1/2,* (*) 



^i+l/2,feW 



(17) 



°j+l/2,k 



(t) — Sj+i^ixj+i/2, Ykj*) 

AX 



S 



j+l/2,k 



Sj+l,k(t) 



(S x ) j+hk (t) and 



(t) = S j>k (x 



Sj,k\t) 



j+1/2, yk? t) 

AX 



(Sx)j,k(t) 



are the corresponding right and left intermediate values of 
S(x,t K ) at (xj + i/ 2 ,yjfc). 

The local speed of wave propagation aJ +1 / 2 * s estimated 
at the cell boundaries (2^+1/2, 2/fc) as the upper bound 



l/2,fc(*) = maX { I^J+l/2,fc (t) /'(<*>) I } 



where a; is a value between S^ +l ^ 2k (t) and Sj + x/ 2 The 
velocity field used in the KT scheme is obtained directly from 
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the Raviart-Thomas space on the cell edges: 

X V j+1 / 2 ,k{t) ■= {Vr) jk {t), ^-l/2,fe(t) := {Vl) jk {t), 

where v r and v\ stand for the velocity on the one-dimensional 
"right" and "left" faces of the cells, 
(ii) Analogously, compute the difference of the one-dimensional 
numerical flux in the y direction keeping x constant and equal 
to Xj. This difference is denoted by 



(19) 



(20) 



^3,^+1/2 



it) :-. 



TTV 



AY 



The one dimensional numerical flux in the y direction is 
1 



H v 

J,fc+l/2 



(t) 



2 L 



V+1/2W f(S+ k+1/2 (t)) 

+ y Vj,k+l/2(t) f{S] 

,k+l/2 (*)) 



a j,fc+l/2( t ) 



si 



>(t) 



J j,k+l/2\ v ) w j,fc+l/2^ 

In a similar way, the correspondent "up" and "down" inter- 
mediate values of S(x,t K ) at (xj,y k +i/2) are 

Ay 



s 



j,k+l/2 



(t) = S jih+1 (t) 



■(Sy) j; k+i(t) and 



^+1/2(0 



- Ay 



The local speed of wave propagation a v - k+1 / 2 (t) in the y di- 
rection is estimated at the cell boundaries (xj, yk+1/2) as t ne 
upper bound 



y 

a j+l/2,k 



:max{|ViAfc(*) /»!} 



where u> is a value between Sf k+1 , 2 (t) and S~ k+1 , 2 (t). Analo- 
gously the velocity field in the y direction is obtained directly 
from the Raviart-Thomas space on the cell edges: 

where v u and Vd stand for the velocity on the "upper" and 
"lower" faces of the cells. 
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(iii) The cell average S • k in the next time step t™ + At c is then 
the solution of the following differential equation 

j^it) = -(^ 1/2it (t) + ^ 1/2 (t)) 

{ ' AX AY 

The numerical derivatives are computed using the MinMod limiter 
given by Equation (fT4l) . In our numerical experiments, the parame- 
ter 9 assumes values 1 < 9 < 1.8. 

The two-dimensional semi-discrete formulation (|2"T|) comprises a 
system of nonlinear ordinary differential equations for the discrete 
unknows {Sj t k{t)}- To solve it, we integrate in time introducing a 
variable time step At 11 . Although the forward Euler scheme can be 
used, it may be advantageous to use higher order discretizations in 
numerical simulations. The numerical examples presented below use 
third-order Runge-Kutta ODE solvers based on convex combinations 
of forward Euler steps. See [12] and [13] for more details on a whole 
family of such schemes. 



4. TWO-DIMENSIONAL NUMERICAL EXPERIMENTS 

We present and compare the results of numerical simulations of 
two-dimensional, two-phase flows associated with two distinct flood- 
ing problems using the KT and NT schemes. 

In all simulations, the reservoir contains initially 79% of oil and 
21% of water. Water is injected at a constant rate of 0.2 pore volumes 
every year. The viscosity of oil and water used are /i = 10.0 cP and 
fi w = 0.05 cP. The relative permeabilities are assumed to be: 

where s ro = 0.15 and s rw = 0.2 are the residual oil and water satu- 
rations, respectively. 

For the heterogeneous reservoir studies we consider a scalar abso- 
lute permeability field K(x.) taken to be log-normal (a fractal field, 
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see [6] and [21] for more details) with moderately large heterogene- 
ity strength. The spatially variable permeability field is defined on a 
256 x 64 grid with three different values of the coefficient of variation 
CV (CV = 0.5, 1.2, 2.2) given by the ratio between the standard 
deviation and the mean value of the permeability field. 

We now discuss the simulations in the slab geometry. We consider 
two-dimensional flows in a rectangular, heterogeneous reservoir (slab 
geometry) having 256m x 64m with three different sizes of computa- 
tional grid: 256 x 64, 512 x 128 and 1024 x 256 cells. The boundary 
conditions and injection and production specifications for the two- 
phase flow equations ([Q) ~ ( I^ZI ) £ir<3 cLS follows. The injection is made 
uniformly along the left edge of the reservoir and the production 
is taken along the right edge; no flow is allowed along the edges 
appearing at the top and bottom of the reservoir. 

Figures [HI HI and [5] refer to a comparative study of the NT and 
KT schemes, showing the water saturation surface plots after 275, 
250 and 225 days of simulation for the three different CV values 
(CV = 0.5, , 1.2, and 2.2). The results obtained with the NT scheme 
were computed using three computational grid: the coarsest grid 
with 256 x 64 cells, and two levels of refinement denoted by NTr 
and NTrr with 512 x 128 and 1024 x 256 cells, respectively (See the 
first three pictures from top to bottom of Figures [3[ HJ and |3j). At 
the same time, the bottom pictures in those figures are the results 
presented by the KT dimension by dimension scheme on the coarsest 
computational grid of 256 x 64 cells. 

For each heterogeneity, we computed the difference between the 
results produced by the NT scheme in the three computational grids 
with the corresponding result produced by the KT scheme in the 
coarsest grid. We consider the solution of the NT scheme in the finer 
grid as the reference solution. The differences are then computed 
using the L 2 norm relative to this reference solution as follows 



II F - G 

;22) numerical differences = — : 

' ; NTrr 
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Here F stands for the KT solution and G for a NT solution. The 
graph in Figure [2] shows these differences. Note that as we refine 
the solution of the NT scheme, the differences become smaller indi- 
cating a comparable accuracy for the simulations produced by the 
KT scheme in the coarsest grid and the result produced by the NT 
scheme in the finest grid. These differences indicate that one has to 
refine twice the grid used in the NT scheme to produce an equivalent 
solution to the one produced by the KT scheme using the coarsest 
grid. 

We now turn to the discussion of the set of simulations performed 
in a five-spot pattern. In case of a five-spot flood discretized by 
a diagonal grid (Figure [6]), injection takes place at one corner and 
production at the diametrically opposite corner; no flow is allowed 
across the entirety of the boundary. In case of a five-spot flood 
discretized by a parallel grid (Figure EJ), injection takes place at two 
opposite corners (say, bottom left and top right), and production is 
through the remaining two corners (say, bottom right and top left). 
Figures [6] (diagonal grid) and [7J (parallel grid) show the saturation 
level curves after 260 days of simulation obtained with the NT and 
KT schemes for two levels of spatial discretization. 

In both Figures [6] and the pictures on the left column are the 
results obtained with the NT scheme and the ones on the right were 
computed with the KT scheme. In these Figures, the grids are re- 
fined from top to bottom and have 64 x 64 and 128 x 128 cells in 
the diagonal pattern and 90 x 90 and 180 x 180 cells in the parallel 
grid. 

It is clear that the KT scheme (right column pictures in Figure [6] 
and in Figue [7J is producing incorrect boundary behavior. More- 
over, as the computational grid is refined (right column and bottom 
picture in Figures [6] and UJ) this problem seems to be emphasized. 
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A CV=0.5 ♦ CV=1.2 o CV=2.28 




0.02 

NT(256x64) NTr(512x128) NTrr( 1 024x256) 

Figure 2. The numerical differences in L 2 norm be- 
tween the solution of the KT scheme and the solutions 
of the NT scheme using three computional grids. As 
we refine the grid of NT scheme, the differences be- 
come smaller. 
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Figure 3. Water saturation surface plots after 275 
days of simulation in a heterogeneous reservoir having 
256 m x 64 m, with CV = 0.5 and viscosity ratio 
20. The first three pictures from top to bottom used 
the NT scheme with grids having 256 x 64, 512 x 128 
and 1024 x 256 cells, repectively. The bottom picture 
shows the KT scheme with a grid of 256 x 64 cells. 
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Figure 4. Water saturation surface plots after 250 
days of simulation in a heterogeneous reservoir having 
256 m x 64 m, with CV = 1.2 and viscosity ratio 
20. The first three pictures from top to bottom used 
the NT scheme with grids having 256 x 64, 512 x 128 
and 1024 x 256 cells, repectively. The bottom picture 
shows the KT scheme with a grid of 256 x 64 cells. 
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Figure 5. Water saturation surface plots after 225 
days of simulation in a heterogeneous reservoir having 
256 m x 64 m, with CV = 2.2 and viscosity ratio 
20. The first three pictures from top to bottom used 
the NT scheme with grids having 256 x 64, 512 x 128 
and 1024 x 256 cells, repectively. The bottom picture 
shows the KT scheme with a grid of 256 x 64 cells. 
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(a) NT: 64x64 grid (b) KT: 64x64 grid 




8 16 24 32 40 48 56 OB 16 24 32 40 4B 56 64 



(c) NT: 128x128 grid (d) KT: 128x128 grid 

Figure 6. Water saturation level curves for two- 
phase flows in a five-spot well configuration - diagonal 
grid. 
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X 



(c) NT: 180x180 grid (d) KT: 180x180 grid 



Figure 7. Water saturation level curves for two- 
phase flows in a five-spot well configuration - parallel 
grid. 
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5. Conclusions 

In one spatial dimension the KT scheme is a small modification of 
the NT scheme which uses more precise information about the local 
speed of propagation. This approach leads to a very simple numer- 
ical recipe producing numerical solutions more accurate then those 
provided by the NT scheme. On the one hand, In two spatial dimen- 
sions the KT scheme uses numerical fluxes in the x and y directions 
that can be viewed as generalizations of one-dimensional numerical 
fluxes. This is called the dimension by dimension approach. The 
NT, in the other hand, uses a genuinely two-dimensional configu- 
ration. In the case of a slab geometry, the fluid flows mostly in 
one direction. For this reason, this flow may be viewed as a one- 
dimensional flow and the KT scheme is expected to produce very 
accurate solutions like those presented in Figures [31 0] and [5j In 
the five-spot problem, the fluid flows in both x— and y— directions, 
causing a genuinely two-dimensional displacement. The KT scheme 
produces incorrect boundary behaviors in the five-spot numerical ex- 
amples. We remark that this incorrect behavior is not present in the 
results produced by the NT scheme. The dimension by dimension 
approach of the KT scheme might be a source of numerical errors 
for a class of problems with an intrinsic two-dimensional geometry. 
These numerical errors may lead to incorrect behavior like those in 
the five-spot problem. The authors are currently working on an im- 
provement of these schemes in order to compute more precisely a 
genuinely two-dimensional numerical flux. 
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